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ABSTRACT 

We examine the distribution of neutral hydrogen in cosmological simulations carried 
out with the new moving- mesh code AREPOand compare it with the corresponding 
GADGET simulations based on the smoothed particle hydrodynamics (SPH) technique. 
The two codes use identical gravity solvers and baryonic physics implementations, but 
very different methods for solving the Euler equations, allowing us to assess how nu- 
merical effects associated with the hydro- solver impact the results of simulations. Here 
we focus on an analysis of the neutral gas, as detected in quasar absorption lines. We 
find that the high column density regime probed by Damped Lyman-a (DLA) and 
Lyman Limit Systems (LLS) exhibits significant differences between the codes. GAD- 
GET produces spurious artefacts in large halos in the form of gaseous clumps, boosting 
the LLS cross-section. Furthermore, it forms halos with denser central baryonic cores 
than AREPO , which leads to a substantially greater DLA cross-section from smaller 
halos. AREPO thus produces a significantly lower cumulative abundance of DLAs, 
which is intriguingly in much closer agreement with observations. For the low column 
density gas probed by the Lyman-a forest, the codes differ only at the level of a few 
percent, suggesting that this regime is quite well described by both methods, a fact 
that is reassuring for the many Lyman-a studies carried out with SPH thus far. While 
the residual differences are smaller than the errors on current Lyman-a forest data, 
we note that this will likely change for future precision experiments. 
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quasars: absorption lines - methods: numerical 



1 INTRODUCTION 

Absorption features in quasar spectra offer a unique view 
of structure formation at redshifts z — 2 — 4. Lyman- 
a absorption of neutral hydrogen directly tracks the dis- 
tribution of gas during the initial stages of galaxy forma- 
tion. Damped Lyman-a Systems (DLAs) ha ve a neutral 



hydrogen column density Nm > 10 cm (Wolfe et al. 



1986), and thus can be observed through natural broad- 



ening of the Lyman-a line. DLAs are thought to be high 
redshift proto-galaxies, sufficiently dense that their interi- 
ors are shielded from the ionising effect of the diffuse ra- 
diation background ( Katz et al.||1996 ). Thus, at z — 2 — 4 



* E-mail: spb@ias.edu 



they are understood to be reservoirs containing most of the 
neutral hydrogen in the Universe 



(e.g. 



et al.||1998| |Gardner et al.||20011 |Maller et al.|20oi| lOkoshi 



ft Nagashima |2005| |Nagamine et al.||2007| |Pontzen et al 



2008| |Barnes ft Haehnelt|2009| |Hong et aI.|20To 



Katz et al 



1996 



Gardner et"aT|[T997l [Prochaska ft Wolfe||1997|" Haehnelt 



Nagamine 



20111 \Cen 



et al | [20101 |McQuinn etaLl[20lT] |Altay et aL 
2012|) . A wide range of observational surveys (e.g. |Wolfe 



et al.||1995| |Storrie-Lombardi ft Wolfe 2000] |Rao ft Turn- 



shek||2000| |Peroux et al.|2001| |Prochaska et al.|2001| |Chen 



ft Lanzettal |2003| |Peroux et al.||2005| |Prochaska eTHI 
2005| |Q'Meara et al ||2007| |Prochaska ft Wolfe||2009| |No- 
terdaeme et al.||2009| ), have studied the properties of DLAs 
together with their lower column density cousins, Lyman 
Limit Systems (LLS). LLS are defined to be absorbers with 
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10 17 cm _2 < N m < W cm and have been connected to 



filamentary structures on the outskirts of proto-galaxies ( Fu- 
magalli et al. 2011 Faucher-Giguere &; Keres|2011 ). At lower 
column densities, we see the Lyman-o; forest; a complex re- 
gion of overlapping Lyman-a lines. The Lyman-oj forest is 
a probe of the matter distribution in the low-density in- 
tergalactic medium (Hernquist et al. 1996; Miralda-Escude 



et aI.|1996||Dave et al.|1999| |Croft et al.|1998| ). It has been 
used to constrain the initial conditions of the Universe, the 
processes that govern the thermal state of the gas, and in- 



directly helium reionisation (e.g. Seljak et al. (2005); Viel & 
Haehnelt|(|2006t;|Bolton et al.|([2008|);|Faucher-Giguere et aT 



( 2008a b ) ;|Viel et al.| 02009) ; |Becker et al.|fl2011| ); |Lidz etaL 
(2010|;|Bird et al.|||2011^ 



The Lyman-a forest is produced by diffuse absorbers 
~ 100/i _1 kpc across, collapsing under gravity. DLAs, on 
the other hand, are produced in smaller, denser systems 
and are strongly influenced by gas physics. In both cases, 
obtaining accurate quantitative results for their properties 
requires following non-linear gravitational collapse with N- 
body simulation techniques and a method for solving the 
inviscid Euler equations for the cosmic gas. Historically two 
main approaches have been used when tackling this prob- 
lem through direct simulations: grid-based codes that dis- 



cretise the gas on an adaptively refined mesh (Berger & 
iColellaligsgllTeyssierpOO^llO'Shea et al.|2004| and particle- 
based codes that use the smoothed particle hydrodynamics 



(SPH) technique (|Lucy||1977| |Gingold fe Monaghan||1977 
Monaghan 1992, 2005). The former solves the equations of 



motion for a fluid on a stationary grid of cells in an Eulerian 
fashion. The latter is a pseudo-Lagrangian technique, where 
the fluid is split into a number of discrete mass elements, 
which are assumed to be indivisible and are then followed 
as particles. Smoothed fluid quantities are constructed from 
the particles through kernel interpolation. 

Each approach has its own advantages and disadvan- 
tages, making it non-trivial to judge their relative accuracy 
for different applications. For example, while SPH codes 
have excellent conservation properties, they suffer from large 
gradient errors and accuracy problems in resolving fluid in- 
stabilities ( |Agertz et al. 2007), something that can make 
them fail dramatically in some idealised fluid dynamics prob- 
lems ( |Sijacki et al.|20lT ). Eulerian mesh codes on the other 
hand offer high accuracy for capturing shocks, but many 
have incorporated relatively inaccurate gravity solvers (see, 
e.g. O'Shea et al.|2 005). Furthermore, their truncation error 
depends on the absolute velocity of the gas, which is prob- 
lematic for the often highly supersonic flows encountered in 
cosmological simulations. 

Springel (2010) proposed a new technique which aims to 
combine the strengths of both approaches. In this method, 
space is tessellated with a moving mesh; cells advect with the 
fluid flow, and as a result each cell contains approximately 
the same gas mass. Thus the resolution automatically in- 
creases in areas of higher density, just as in Lagrangian 
methods like SPH. A moving mesh method can resolve 
shocks equally as well as grid methods, since in both cases 
the Euler equations are solved with a high-accuracy finite 
volume method, but without the advection errors present in 
Eulerian codes. The AREPO code is an implementation of 
this technique and has been shown to perform well in many 



idealised simulations of fluid dynamical problems (Springel 
2010} jSijacki et al.|2011 



We compare results obtained with AREPO to corre- 
sponding simulations with the well-tested SPH code GAD- 
GET , last described in |Springel| ( |2005| . Although we know 
that SPH cannot accurately describe fluid instabilities and 
mixing processes in its widely employed 'standard' form, it 
is not obvious a priori to what extent these inaccuracies 
are important for different aspects of structure formation, 
or how they affect our interpretation of observational data. 
This paper therefore examines the important issue of how 
predictions for the neutral hydrogen distribution in the uni- 
verse depend on the employed numerical technique. 

Our study is part of a series which compares simulations 
run using AREPO and GADGET with identical initial condi- 
tions, gravity solver, and baryonic physics parameters. They 
thus differ only in the approach to hydrodynamics, offering 
an unrivalled opportunity to isolate the effects of numerical 



uncertainties. In previous work, Vogelsberger et al. (2011) 



examined the global properties of baryons and halos, as well 
as performing convergence and numerical tests on AREPO . 



Sijacki et al. (2011 ) studied a number of idealised fluid prob- 



lems to clarify the origins of the observed differences. |Keres 



et al. (2011) looked at the effect on the resulting galaxy 
properties, and Torrey et al.| ( |2011| focussed on the struc- 
ture of galactic discs. Here we shall look at how the use of a 
moving- mesh technique affects the properties of DLAs, LLS 
and the Lyman-a forest. Neutral hydrogen is particularly 
relevant as it is a comparatively clean probe of the hydrody- 
namics. We shall include star formation, but neglect strong 
feedback from outflows. Also, rather than a full radiative 
transfer model, we shall use a simple density cut-off to ac- 
count for self-shielding of the gas. While less sophisticated 
than many previous works, this simple approach allows us to 
focus on the effect of the hydro solver and helps to connect 
our intuition from idealised tests with observations. 

We build on a large literature of simulated DLAs: 
Nagamine et al. ( 2004a|b ) used SPH simulations with GAD- 
GET to examine their abundance and metallicity. They also 
examined the effects of galactic winds on DLAs, looked at 
further by|Nagamine et al.|(|2007|) andjTescari et al.|([2009). 



Pontzen et al. ( 2008 ) attempted to reproduce many observed 



properties of DLAs with a complete simulation framework 
incorporating a simple model of radiative transfer and su- 
pernova feedback. Radiative transfer was also examined in 



more detail by jYajima et al.| fl2011| ). |Cen| ( |2012| ) and [Fu 



magalli et al. (2011) both used Eulerian grid codes with 



adaptive mesh refinement to study DLAs and LLS. They 
included feedback from star formation and metal cooling, 
which we neglect, but their AMR-based simulations were 
unable to resolve halos less massive than 10 10 h~ 1 M(7 ) . 

This paper is structured as follows. In Section [2] we 
discuss our methods in more detail. We then present our 
results in Section [31 and conclude with a summary of our 
findings in Section^] 



2 METHODS 

2.1 Numerical codes and simulation set 

Both AREPO and GADGET compute gravitational interac- 
tions using the TreePM approach, as described in Springel 
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(2005). Radiative cooling is implemented, following Katz 
et al.| ( p~996| , using a rate network, including line cooling, 
free-free emission and inverse Compton cooling off the cos- 
mic microwave background. We assume ionisation equilib- 
rium and optically thin gas during the simulation, account- 
ing for gas self-shielding by post-processing the simulation 
outputs, as described in Section |2?2| The history of the UV 
background follows the estimates of Faucher-Giguere et al.| 
([2009}. 

Star formation is implemented with the effective two- 
phase model of Sprin geT &; Hernquist| ([2003 ) . We use the 
same parameters as in that paper, giving a threshold density 
for star formation of Nh = 0.13cm -3 . The star formation 
time-scale is assumed to scale with density as £* oc p -0-5 , 
normalized to 2.1 Gyr at the threshold density in order to 
match the local relation between the gas surface density and 
star formation rate of galaxies ( Kennicutt|1998| ). The energy 
associated with supernovae heats the multi-phase medium 
and thereby regulates star formation, but in the model im- 
plementation used here this feedback is not strong enough 
to drive significant outflows. We note that McDo nald et al.| 
( 2005 ); Cen et al. ( 2005 ) found that galactic winds have little 



effect on the Lyman-a forest, but Nagamine et al. (2007); 
Tescari et al. (2009) showed that they do affect the DLA 



cross-section and column density function. Since this work is 
primarily aiming to compare the relative performance of two 
different codes, the absence of strong winds is not a prob- 
lem. Indeed, it helps to avoid more complicated gas motions 
which could obscure our results. 



We use the simulations first described in Vogelsberger 
|et al.| ( [20TT| . The initial conditions were generated at z — 99, 
using the power spectrum fit of |Eisenstein &; Hu| ([1999 ) with 
cosmological parameters Q m = 0.27, Qa = 0.73, = 0.045, 
cr 8 0.8, n s 0.95 and H 70kms~ 1 Mpc~ 1 (h 0.7), 



consi stent with the latest WMAP results (Komatsu & et al. 
2011). The box size is 20/i _1 Mpc. 



Most of our results are from the highest resolution sim- 
Vogelsberger et al.| ( |201l| ). This has 512 3 dark 

'Mo, 



ulation of 



matter particles, with a particle mass of 3.722 x 10 h~ 
and a comoving gravitational softening length of 1 /i _1 kpc 
(1/40 of the mean inter-particle spacing). Each simulation 
was initialised with 512 3 gas elements, with an initial mass 
of 7.44 x 10 5 /i _1 M©. However, the number and mass of gas 
elements are both altered over time by star formation, and, 
in AREPO , mass fluxes, refinement and de-refinement of 
grid cells. Further details on the refinement implementation 
used in AREPO may be found in Vogelsberger et al. (2011 ). 



We examined three output times corresponding to redshifts 
z = 4, 3 and 2. We emphasize that all the simulations used 
the same realisation for their initial conditions as well as the 
same parameters for the sub-grid physics model, allowing a 
comparison on an object-by-object basis. 



We refer the reader to Springel (2010) for further details 



of the moving mesh implemented in AREPO , to Springel 



(2005) for details of the gravity computation, to Springel & 



Hernquist|(|2002| for the SPH imp lementation in GADGET , 
and to | Vogelsber ger et al.| ( |2011| for a full account of the 
parameters of the simulations. In the rest of this Section, we 
focus on those aspects of the analysis specifically concerned 
with DLAs. 



2.2 Gas self-shielding 

The gas in DLAs is self-shielded, with a neutral fraction 
close to unity. A full treatment requires radiative transfer 
and has been studied in, e.g., Pontzen et al. (2008); AL 



tay et al. (2011); Yajima et al. (2011); Cen (2012). Yajima 



et al. (2011), using outputs from an SPH simulation post- 



processed with radiative transfer, suggested modelling the 
transition to self-shielded gas with a step function. We con- 
sidered this, but found that it caused an artificial kink in the 
column density function, especially prominent for AREPO , 
although it did not affect any of our other results. Instead we 
fit the neutral fraction as a function of density to the results 
of their preferred simulation. Gas is assumed to transition 
between equilibrium with the UVB for p < pu and complete 
self-shielding for p > p s , with the transition region given by 



n^Y B (Ps-p) p + (p-pur 

(ps - p\j) P 



(1) 



where n^i 75 is the neutral fraction when in equilibrium with 
the UVB. The best fit values were p = 2.68, p v = 4.53 x 
10~ 3 cm -3 and p s — 1.52 x 10~ 2 cm -3 . Furthermore, we 
verified that this was consistent with the results of Altay 
|et al.| (2011 ), and that the DLA abundance is not sensitive 
to the exact position of the transition. 



Nagamine et al. (2004a) proposed identifying self- 
shielding with the onset of star- format ion; i.e., p ss = 
0.1289 cm -3 . This is an order of magnitude higher than the 
value of p ss suggested by radiative transfer and produces a 
slight, probably unphysical, coupling of the neutral fraction 
to the star formation rate. While we did not use this pre- 
scription, we did check its effect and found that, although 
the absolute value of many DLA properties changed (which 
is to be expected when changing the density threshold for 
self-shielding by an order of magnitude), the differences be- 
tween AREPO and GADGET were largely unaffected. This 
gives us further confidence that our results are robust to 
changes in the self-shielding prescription. 



Cen (2012) and Altay et al. (2011) included a prescrip- 



tion for the formation of molecular hydrogen. | Altay et al.| 
( |20TTt were motivated by |Blitz &; Rosolowsky| ( |2006| , who 
found that the surface density of molecular hydrogen in local 
galaxies was strongly correlated with the hydrostatic pres- 
sure. Altay et al. (2011 ) assumed the same relation held be- 



tween molecular hydrogen density and pressure in the star- 
forming phase of the ISM at high redshift. We considered 
such a prescription, but found the only noticeable effect was 
on the column density function for Nm > 10 22 atoms cm -2 , 
and was not sufficient to bring it into agreement with obser- 
vations. Furthermore, we compared the shape and amplitude 
of the molecular hydrogen column density function predicted 
by the model at z = to the observed values of Zwa an fc] 



Prochaska (2006) and found that they did not match. We 



believe both these facts are because the lack of feedback in 
our simulations led to a surplus of over-dense material. We 
therefore decided not to incorporate molecular hydrogen in 
our analysis until our simulations include strong feedback 
processes. 



© 2011 RAS, MNRAS 000,[T}[T3] 



4 S. Bird et al 



2.3 DLA selection and column density 

We use a halo catalogue generated with the Friends-of- 
friends (FOF) algorithm ( |Davis et al.||1985| and a linking 
length of 0.2 of the mean inter-particle spacing. The FOF 
algorithm is applied only to the DM particles, whereas bary- 
onic particles/cells are later assigned to their nearest DM 
particle and included in the halo to which the correspond- 
ing DM particle belongs. Self-bound concentrations of mass 
are identified within each FOF halo using the SUBFIND al- 
gorithm (Spri ngel et al.|2001 ), modified to account correctly 
for baryons ( |Dolag et al.|2009| . SUBFIND identifies halo sub- 
structure by generating an adaptively smoothed density field 
and searching for gravitationally bound over-densities. Fol- 
lowing a common convention in the literature, we use M 2 oo 
as an estimate of the halo virial mass. M200 is the mass en- 
closed within a spherical region of radius R200 within which 
the mean density is 200 times the critical density. 

We consider only resolved halos with M > 
400(Q m /Qb) rnb, where vrtb is both the SPH gas particle 
mass and the initial mass of AREPO cells. For our 2 x 512 3 
simulation, this implies a halo resolution threshold of M ^ 
2 x 10 9 hT 1 Mq. We note that our conclusions are unaffected 
by the exact placing of this limit. To avoid confusing over- 
laps, we remove halos with larger neighbours closer than 
R200 from our catalogue. 

To find the projected neutral hydrogen column density 
around each halo, we consider a grid of size 2i?2oo cen- 
tred on each halo, divided into equal-sized cells with side 
length equal to the gravitational softening length (for us 
~ l/i _1 kpc). Gas elements in the halo are projected onto 
this grid, using an SPH kernel. In the case of AREPO , the 
SPH smoothing length is chosen so that the volume covered 
by the SPH kernel is identical to the volume of the mov- 
ing mesh cell. We also considered a cloud-in-cell kernel and, 
like |Nagamine et al.] ( |2004a| , found that it made negligible 
difference to our results. 

We calculate the column density along a line of sight 
by projecting our interpolated density field along a single 
direction, here the x axis. The column density is given by 



iVHi 



S^e(l + ,) 2 ; 



nip 



(2) 



where mp is the proton mass, e is the side length of a single 
grid cell and pm(x) is the average neutral hydrogen density 
inside the cell. The factor of (1 -\- z) 2 enters because pm is in 
comoving units and Nm is in physical units. We checked dif- 
ferent projection directions and found the overall statistical 
properties of DLAs were unchanged. However, since galaxies 
in AREPO are more disc-shaped, they tended to look some- 
what different when viewed edge-on. We also verified that 
our results are robust to changes in the resolution of the grid 
used for the map making. 



2.4 Column density function 

The column density function, f(Nm), is defined observa- 
tionally such that f(Nm) dNm dX is the number of ab- 
sorbers per sight-line with column density in the interval 
[iVm, iVm + diVHi]. We identify sight-lines with grid cells and 
thus count absorbers by computing a histogram of the col- 
umn density on the grid. Equating grid cells with sight-lines 



assumes that two simulated DLAs will rarely be found along 
the same sight-line, which is a good assumption given the 
small size of our box. More explicitly, we define the column 
density function by 



f(N) 



F(N) 
AN 



AX(z) . 



(3) 



F(N) is the fraction of the total number of grid cells in 
a given column density bin, and AX(z) is the absorption 
distance per sight-line. As described by Bahcall & Peebles 
(1969); Naga mine et al.| ( |2004a| , the (dimensionless) absorp- 
tion distance is given by 



(4) 



x ^ = f (1 + z ' )2 *M d *'' 

and, for a box of comoving length AL, we have AX 
(H /c)(l + z) 2 AL. 



3 RESULTS 

In this section, we present our results for the comparison 
between AREPO and GADGET . In Sections [3j] and we 
look at the effects on large and small halos, respectively. 
Then we examine various statistical properties of the DLAs: 
Sect ion 1 3 .3| considers the DLA and LLS cross-sections, Sec- 
tion 13.41 the observed DLA abundance and Section 13.51 the 
column density function. Finally, in Section |3.6| we discuss 
our results for the Lyman-a forest. We emphasize that we 
have checked explicitly that our results are unchanged when 
using simulations with 2 x 256 3 particles - a factor of 8 lower 
mass resolution - at least for halos above the resolution limit 
of the lower resolution simulations. 



3.1 Large halos 

Figure [l] shows the distribution of neutral hydrogen in the 
largest halo in our simulation for z = 4 — 2. For column 
densities with A^hi > 10 19 atoms cm -2 , gas is almost en- 
tirely self-shielded and thus is traced extremely well by the 
neutral hydrogen. We can see that GADGET produces a 
large number of small, circular, gaseous artefacts, which are 
largely absent in AREPO . Similar "blobs" have been seen 
in other cosmological SPH simulations (including calcula- 
tions with different SPH codes), and we interpret their exis- 
tence as a numerical artefact due to the suppression of fluid 



instabilities and mixing in SPH; see |Torrey et alT (2011) 
for further details. In GADGET, these SPH blobs make 
up the bulk of the LLS, especially in the column density 
range 10 17 < N m < 10 19 atoms cm -2 . AREPO reveals that 
LLS are more commonly produced from distinct filamen- 
tary structures. Figure [2] shows this more explicitly. Here 
we have discretised the column density map; any cell with 
Nm < 10 17 cm -2 is shown in blue, cells producing LLS are 
yellow, while DLAs are red. 

There is a more subtle change in the DLA cross-section. 
DLAs in AREPO are more concentrated in the centre of 
the halo, but the overall cross-section is not significantly 
changed; although there are fewer substructures, this is par- 
tially compensated by the higher accretion rate of the central 
halo. 
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Figure 1. Neutral hydrogen distribution around the largest halo, at three different redshifts. Each row shows a different redshift, from 
top to bottom, z = 4, 3, and 2. The left column shows AREPO and the right GADGET . y and z are comoving coordinates. 



3.2 Small halos 

Figure [3] shows the radial profile obtained by stacking all 
halos with mass 3 x 10 9 /i _1 M < M < 3.5 x 1O 9 /i _1 M 
at z = 3. We radially integrated the neutral hydrogen grids 
shown in Figure [I] and averaged over all halos in the given 
mass range. These halos are approximately spherical, so the 
effect of substructure is relatively small. The column density 
in the central 5/i _1 kpc of the halo is larger by a factor of 
about 7 in GADGET . This effect is more pronounced at 
redshift z — 4 (a factor of 10) than at z — 2 (a factor of 3). 



Its origin does not appear to be numerical; the density of the 
innermost 5/i _1 kpc is larger in GADGET for all halos with 
mass up to M = 5 x 10 10 /i _1 M , in both the 2 x 512 3 and 
2 x 256 3 particle simulations. However, for M > 10 10 /i _1 M 
(and at z = 2), the characteristic size of a DLA is much 
larger than 5 /i _1 kpc, so changes in the central density have 
a much smaller effect. 
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Figure 2. Distribution of gas around the most massive halo at z = 4, highlighted so that all DLA cells are in red and all LLS cells are 
in yellow. Lower density cells are in blue. The left panel shows AREPO , the right panel GADGET . y and z are again comoving. 
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Figure 3. Stacked radial density profiles around the centre of 
halos with mass 3 x 10 9 /i _1 M© < M < 3.5 x 10 9 /i _1 M© at 
redshift z = 3. The blue dashed lines show 662 stacked halos 
from GADGET , while the red solid lines show 719 stacked halos 
from AREPO . The black dot-dashed line indicates the density 
required for being over the DLA cut-off, A^hi = 10 20-3 cm - 2 . R 
is comoving, but the radial density is in physical units. 
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2 


Gadget 


0.849 


31.5 


-50.0 
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Table 1. Numerical parameters of the fit described in Eq. j5}. 



3.3 DLA and LLS cross-section 

We define the DLA cross-section, ctdla, of a halo to be the 
area covered by all grid cells with column density Nm > 



1(T 



atoms cm . The LLS cross-section, ctlls, is defined 



Figure [4] shows the DLA cross-section, ctdla, as a func- 
tion of halo mass. The shaded region delineates the area in 
the M — a dla plane populated by halos. Symbols with error 
bars show the median, upper and lower quart iles of ctdla 
for halos in seven halo mass bins. In order to fit the features 
at small and large halo masses, we modelled ctdla and ctlls 
with a three-component power law 



M 
Mo 



M 



10^/ 5 \M 



10 (0-N)/5_ 



(5) 



where M = 10 10 5 /i _1 M and N = 20.3 for DLAs or 
N = 17 for LLS. The free parameters a, b, c, and e are 
foun d by a simultaneo us fit to ctdla and ctlls, using MP- 
FIT dMarkwardt | 2009f ] The resulting numerical values are 
listed in Table ^ and the fit shown by the lines in Figure 
4) Th is procedure is similar to that outlined in Nagamine 
et al.| ( |2004a| , except that we have excluded halos which do 
not form DLAs from our fit, instead of assigning them an 
arbitrary small cross-section as in that work. 

Qualitatively, our results for ct D la in GADGET are 
in good agreement with those of Nagamine et al. (2004a) 
for overlapping halo masses. They found that, for SPH, 
a turnover in the dDLA-halo mass relation occurred for 
M ~ 10 8-5 /i -1 M0, below the resolution limit of our simula- 
tions. They also found a slight excess in ctdla over their bare 
power law for large halos, which we have fit for explicitly. 
Our ctdla has a somewhat smaller amplitude than the 03 



model of Nagamine et al. (2004a) and the aDLA-halo mass 



similarly, but with 10 20 3 > N m > 10 atoms cm 



relation is somewhat shallower in slope. This is in agreement 
with the results of |Tescari et al.| ( |2009| ) and is probably due 
to a difference in cosmological parameters; |Nagamine et al.] 
(2004a) used cr 8 = 0.9, while we have erg — 0.8. 

Halos in the mass range 10 11 hr^M® > M > 
10 10 H~ 1 Mq show similar cross-sections in both AREPO and 
GADGET . The extra substructure discussed in Section 13.11 
leads to much larger ctlls in GADGET for halos of mass 
larger than 10 11 h~ 1 Mo, especially at low redshift, because 
of the increased number of massive halos. This effect is 
smaller for ctdla, as discussed above, but still present. 



As ported to python by Mark Rivers and Sergey Koposov 
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DLA cross-section at z=4.0 
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DLA cross-section at z=2.0 
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Figure 4. Left panels: The comoving DLA cross-section, <tdla- Right panels: The LLS cross-section, <t lls . Each row shows a different 
redshift, z = 4 (top), z = 3 (middle), and z = 2 (bottom). Regions containing at least one GADGET halo are shown in blue, while 
regions containing at least one AREPO halo are shown in red. The symbols with error bars show the median DLA cross-section in seven 
evenly-spaced mass bins; the error bars show the upper and lower quartiles. Red triangles are for AREPO , while blue squares are for 
GADGET. The symbols for GADGET have been offset horizontally by 10% for visibility. The lines denote the results of our fit; red 
solid lines are for AREPO , while blue dashed lines are for GADGET . 



Due to the effect discussed in Section [372j ctdla is sig- 
nificantly reduced in AREPO for small halos, producing a 
turnover in the power law relation between mass and DLA 
cross-section. This turnover occurs because halos below a 
certain size have a central column density which is on aver- 



age less than the DLA cut-off. It also occurs in GADGET , 
but at halo masses below the resolution limit of our simula- 
tions ( |Tescari et al.|2009 ). Notice that this reduction in the 
cross-section of low-mass DLAs is not carried over to the 
cross-section of low-mass LLS; if anything, small halos in 
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Figure 5. Cumulative abundance of DLAs per unit redshift as a 
function of total halo mass at z = 3 for AREPO (red) and GAD- 
GET (blue) . The yellow bar shows the observed total cumulative 
DLA abundance of |Prochaska et al.| (2005). 



GADGET have slightly lower <jlls than their counterparts 
in AREPO . This suggests that these small halos have ap- 
proximately the same amount of neutral hydrogen in both 
codes, but that this gas is more concentrated in the central 
peak in GADGET . 

While the overall amplitude of ctdla remains roughly 
constant with redshift, that of ctlls decreases fairly substan- 
tially with time, in both codes. Figure [I] shows the reason 
for this; by z — 2, most of the filaments and streams that 
dominate the LLS cross-section have been swept into halos. 



3.4 Incidence rate of DLA systems 



Following Nagamine et al. (2004aJ), we calculate the inci- 



dence of DLAs by convolving the halo mass function with 
our fit to (7 dla • This allows us to account for dark matter 
halos of smaller mass than the resolution limit of the simu- 
lation. The cumulative number of DLAs per unit redshift is 
then denned as 



dz 



(> M,z) 



dr 
dz 



f 

J M 



dn h 
dM 



ODLAdM , 



(6) 



where dnh / dM is the Sheth-Tormen dark matter halo mass 
function, <7dla is given by Eq. ([5]), truncated at M = 
1O 12 - 5 /i _1 M0 and 

dr _ c 
dz ~ JTo 

defines the line element. 

Figure |5] shows diVoLA/ dz for the AREPO and GAD- 
GET simulations. The yellow band marks an observational 
estimate of the total DLA abundance, 



V / ^m(l + ^) 3 +^ 



, fdN^ A \ 



-0.6 ±0.1. 



(8) 



recovered by Nagamine et al. (2007) from the data of 
Prochaska et al.M2005l. Our results for the GADGET simu- 



lation are similar to the weak wind scenario of Tescari et al.l 



(2009) and produce more DLAs than are observed. They 



showed that this discrepancy could be reduced with the ad- 
dition of feedback processes. We find that the lower DLA 



cross-section of small and large mass halos in AREPO com- 
bine to reduce the cumulative DLA abundance by a factor of 
two, bringing it almost into agreement with the upper limit 
from observations. 



3.5 Column density function 

Figure [6] shows the neutral hydrogen column density func- 
tion, denned in Section |2.4| Our results suggest a signifi- 
cantly shallower column density function than preferred by 
observations. This was also seen by |Nagamine et a l. (2004a) 



andlTescari et al. (20091), who found that adding feedback to 



the simulations produced better agreement. We will address 
this question in future work based on forthcoming AREPO 
simulations with stronger wind feedback. 

Overall, the differences in the column density function 
between the codes are fairly limited; a change of ~ 30% of 
the total. At first glance this may seem puzzling; why have 
the generally smaller AREPO halos made so little difference 
to the column density function? Figure [7] shows the fraction 
of the column density function contributed by halos in three 
mass ranges, corresponding to those in Sections |3 . 1 1 and |3 . 2| 
The column density function has been calculated only for 
those halos with the desired mass and then normalised us- 
ing the total GADGET column density function. This allows 
us to clearly see the differences between the codes. The re- 
duced substructure in M > 10 11 h~ x Mq halos does lead to 
a reduced column density for these halos at Nm ~ 10 20 
atoms cm -2 , but their rarity means they make up only for 
a small fraction of the total column density function. The 
largest change is produced by the reduction in cross-section 
of M < 10 10 h~ 1 M Q halos, although this too is diluted be- 
cause lO lo /i _1 M < M < 10 11 h^M^ halos, which make 
up almost as great a fraction of the total column density, 
are mostly unchanged. This explains why differences become 
larger with increasing redshift; we found in Section [3. 2| that 
the effect on M < 10 10 h~ 1 Mo halos became larger at early 
times. 

Both codes produce a slight kink at 10 19 — 10 20 atoms 
cm -2 . This feature is more prominent in AREPO, and is 
a consequence of the transition between self-shielding and 
UVB ionization equilibrium, made more prominent by the 
smoother gas distribution in AREPO . 



3.6 The Lyman-a forest 

For an analysis of the Lyman-a forest, we generated 16000 
simulated Lyman-a spectra with 1024 pixels each from our 
simulations. The positions of the spectra were chosen at ran- 
dom, and particles/cells were interpolated to the lines of 
sight using an SPH kernel. We verified that using cloud- 
in-cell interpolation for AREPO did not affect our results. 
The optical depth from the absorption due to each par- 
ticle was calculated as described in detail in IBird et al.l 
(2011). In order to avoid contaminating the spectra with 



strong absorbers, we did not apply a self-shielding correc- 
tion here. iBolton et al. (2008) found evidence for an in- 



verted temperature-density relation in the Lyman-a forest, 
so that lower density regions have a higher temperature. 
Mechanisms proposed for reproducing this include helium 



reionisation (McQuinn et al. 2009) or volumetric heating 
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Figure 6. The column density function extracted from GADGET simulations (blue dashed lines) and AREPO simulations (red solid 
lines). Green triangles show the constraints of Noterda eme et al~] ( |2009| at z ~ 3, black squares those of lO'Meara et al.| (2007) at z ~ 3.1, 
and the grey region those of |Prochaska et al.| ( |2010| at z = 3.7. 




JV HI (cm- 2 ) ' ' Redshift 



Figure 7. The contribution to the column density function of 
halos with mass M > 10 11 /i _1 M© (solid), 10 11 h' 1 M© > M > 
10 10 /i^M© (dashed) and M < 10 10 h^M® (dotted), nor- 
malised by the total f(N) for GADGET at z = 3. Blue lines 
denote GADGET and red lines AREPO . The solid black line 
shows the total column density function for AREPO divided by 
that for GADGET . 



Figure 8. The effective optical depth, r e ff , as a function of red- 
shift. Solid red lines show AREPO with 2 x 512 3 particles/cells, 
dashed blue lines show GADGET with 2 x 512 3 particles. Dot- 
dashed blue gives GADGET with 2 x 256 3 particles and dashed 
red shows AREPO with 2 x 256 3 resolution. The bottom panel 
shows the percentage difference relative to the GADGET run 
with 2 x 512 3 particles. 



from blazars ( Puchwein et al.||2012| . As our purpose in this 
paper is a code comparison, we did not attempt to model 
this in our simulations. Our temperature-density relation is 
thus somewhat steeper than observed, with 7^5/3 and 
Tocp 7 . 

We define the transmitted flux as T 



exp(-r) ( |Mc-| 

|Donald et aT||2005| |Viel et al.||2009| ), where r is the opti- 
cal depth. A completely transparent Universe has T — 1. 
Observations have determined the effective optical depth, 
7" e ff = T, where T is the mean transmitted flux, the one- 
dimensional flux PDF (a histogram of the flux from each 
spectral pixel) and the flux power spectrum. We extracted 
all three of these quantities from our simulations and com- 
pared the results of AREPO and GADGET . We checked con- 
vergence using the simulations with 256 3 particles. Figure [8] 
shows the effect of a changing particle number on r e ff . GAD- 
GET converged at the 2% level for r e ff and AREPO to - 4%, 
but this convergence becomes significantly poorer for z > 3.5 



in both codes. Bolton & Becker ( 2009| ) found that this oc- 
curs because at high redshift the Lyman-a transmission is 



dominated by progressively less dense regions, which are less 
well resolved. 

In more detail, the effective optical depth is slightly 
lower in AREPO , by 4% at z = 4, 2% at z = 3 and 1% at 
z = 2. Note that r e s is known observationally to ~ 4% at 
z = 3 (Viel et al. 2009| . This small difference is due to a 
slight increase in the mean temperature, To, of the Lyman- 
a absorbers, denned to be gas with T < 10 5 K and p < p Cl 
where p c is the critical density. T in AREPO is ~ 240 K (2%) 
higher than in GADGET at z = 2 — 4. The Lyman-a forest 
is assumed to be in ionisation equilibrium with the UVB, so 
a higher temperature produces a greater ionisation fraction, 
thus decreasing the effective optical depth. These changes 
are somewhat larger for the lower resolution simulations and 
thus are likely to be reduced further with higher numerical 
resolution. It is the standard procedure, when comparing to 
data, to rescale simulated spectra to have the same mean 
flux as is observed. As we are performing a code compari- 
son, we do not rescale our spectra for the presented results, 
but we checked that it did not significantly affect our con- 
clusions. 
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Figure 11. Two-dimensional mass- weighted histograms of the total gas mass as a function of temperature and density, at z = 3, with 
2 x 512 3 resolution. The colour scale shows the total mass of gas elements found in a given temperature-density bin. 
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Figure 9. The flux PDF for AREPO (solid lines) and GADGET 
(dashed lines), at redshifts z = 3.0 (blue), z = 2.6 (green) and 
z = 2.0 (red). Zero flux corresponds to total absorption. Lower 
panel shows the relative change in percent, with positive numbers 
corresponding to a larger PDF in GADGET . 
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Figure 10. The flux power spectrum for AREPO (solid lines) 
and GADGET (dashed lines), at redshifts z = 3.4 (black), z = 3.0 
(blue), z = 2.6 (green) and z = 2.0 (red). The lower panel shows 
the percentage change, with positive numbers corresponding to 
more power in GADGET . 



Figure [9] shows the flux PDF extracted from our simu- 
lations in three redshift bins, chosen to match observations. 
There are small differences between the codes; for T < 0.6 
the flux PDF is larger in AREPO , but for higher transmis- 
sion regions at z = 3, the trend is reversed and AREPO pro- 
duces a lower flux PDF. Figure ^] shows a similarly sized 
effect on the power spectrum of the flux; AREPO predicts a 
slightly lower power spectrum. For comparison, the typical 
uncertainty in the observed flux PDF is 5 — 7% and 5 — 10% 
in the flux power spectrum. These are of the same order as 
the differences we find here between the numerical codes. 
We attribute these changes in the flux PDF and flux power 
primarily to subtle shifts in the temperature-density distri- 
bution. Figure [TT] shows a mass- weighted histogram of the 
temperature and density of gas elements. We can see that 
although the particles follow a similar temperature-density 



relation overall, AREPO produces a wider spread in temper- 
atures for gas near the critical density; this trend continues 
for higher density gas until p ~ 10/? c - 

This difference was somewhat larger for lower-resolution 
simulations with 256 3 particles, which might naively suggest 
that AREPO is converging more slowly in under-dense re- 
gions than GADGET . However, it is not completely clear 
that the two codes are converging to the same result, be- 
cause the differences we see could well be the effect of weak 
structure formation shocks. While these are followed accu- 
rately in AREPO , they are largely lost in GADGET , due 



to SPH's poorer shock-resolving capability (see, e.g. Keshet 
et al. 20031). We would thus expect that a fully resolved 



AREPO simulation would still produce more gas elements 
scattered off the mean temperature density relation than 
GADGET , consistent with our results. 
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We recall that we have not included feedback in our 
simulations; this is expected to affect the Lyman-a flux PDF 



at the 5% level (Bolton et aLl|2008| |Viel et al.||2012| ) at 
z = 2 — 2.5. Feedback processes such as galactic winds may 
interact with the gas dynamics differently in AREPO than in 
GADGET , potentially adding further systematic differences 
between the codes of a similar magnitude. 



4 CONCLUSIONS 

We compared the statistics of neutral hydrogen absorption 
in the SPH code GADGET with that in the moving mesh 
code AREPO . Our aim was to understand how known differ- 
ences in the treatment of the fluid equations manifest them- 
selves as changes in the observable properties of DLAs, LLS 
and the Lyman-a forest. There were significant qualitative 
differences: in GADGET , DLAs and LLS in large halos were 
primarily produced by small spherical objects, and the only 
difference between the two classes of absorbers was the cen- 
tral density. However, in the AREPO simulations, DLAs and 
LLS came from quite different sources. The bulk of the DLA 
cross-section in a large halo was from a central disc, while 
LLS were produced in filamentary structures. This made lit- 
tle quantitative difference to the DLA cross-section, but it 
suggests a different interpretation of high column density 
systems and could potentially be reflected kinematically in 
detailed line profile shapes. We would suggest that any fu- 
ture studies sensitive to the detailed distribution of neutral 
hydrogen inside massive halos, especially LLS, should avoid 
using SPH. 

We found that GADGET produced clouds of clumpy 
substructure around halos with M > 10 11 h~ x M , which 
were essentially absent in AREPO . These GADGET clumps 
had a peak column density of 10 18 — 10 19 cm -2 and so 
boosted the LLS cross-section of these halos significantly. 
There was a similar, but somewhat smaller, effect on the 
DLA cross-section. Furthermore, halos in AREPO had cen- 
tral densities which were lower than their counterparts in 
GADGET . This significantly lowered the DLA cross-section 
in AREPO for objects sensitive to the density in the cen- 
tral 5/i -1 kpc (in practice, halos with M < 10 10 H~ 1 Mq at 
z > 2), but did not greatly affect the LLS. 

These systematic changes made the halo mass - cross- 
section relation shallower where it was dominated by the 
changes to large halos and steeper where it was dominated 
by changes to small halos. The former occurred for LLS, and 
DLAs at redshift z = 2, and the latter for DLAs at z = 3 
and 4. Furthermore, both changes acted to reduce the over- 
all abundance of DLAs in AREPO . The DLA abundance 
for GADGET simulations is somewhat in excess of the ob- 
served value. This discrepancy can be removed with the in- 
corporation of galactic winds into the simulation ( Nagamine 
et al.||2004a| |Tescari et aT||2009) , but we found that our 
AREPO simulations already substantially reduce it, even 
without strong feedback. This was mostly due to the reduced 
cross-section of small halos; large halos are sufficiently rare 
that changes in their cross-section do not have a large effect 
on the total abundance. It also shows that a calibration of 
feedback strength from the DLA abundance would be com- 
promised by numerical effects in SPH. 

The column density function shows a similar pattern; a 



modest reduction in amplitude for Nm = 10 20 — 10 22 cm -2 
and z > 2, driven by the lower central density of small ha- 
los. We found that the amplitude of the column density 
function for Nm > 10 21 cm -2 is still much larger than 
observed, which can be interpreted as evidence that some 
feedback in the form of outflows is needed. Although the 
column density function for large halos was significantly re- 
duced in AREPO by the changes in their substructure, the 
relative rarity of these halos meant that they did not have 
a strong impact on the total column density function. Note, 
however, that feedback processes reduce the amplitude of 
the high-end column density function by removing preferen- 
tially the high-column density cells in small halos; the more 
massive the halo, the less it is affected by winds (|Tescari 



et al.|2009| . Thus, although it appears as if the reduced sub- 
structure around large halos only produces subtle changes 
in this statistic, once galactic winds have blown the small 
halos away, what remains may be very greatly affected. 

Finally, we examined basic Lyman-a forest statistics 
and found that differences between the codes were typically 
at the few percent level. This was due partly to a slightly 
changed thermal history (which is typically marginalised out 
in cosmological studies of the Lyman-a forest) and partly 
to an increased width in the temperature-density relation of 
the gas in AREPO . These differences are small when com- 
pared to changes in the high column density systems. This 
is not unexpected; the evolution of the Lyman-a forest is 
dominated by gravitational effects, hence issues such as the 
accuracy problems of SPH for fluid instabilities do not play 



an important role (a similar result was found by Regan et al. 
2007). They are however comparable to the statistical un- 



certainties of current data, and thus may bias derived pa- 
rameters. Although a full cosmological analysis is beyond 
the scope of this paper, it seems likely that the largest dif- 
ference will occur in the derived thermal history of the IGM, 
since the change in the flux power spectrum came predom- 
inantly from an alteration in the temperature-density rela- 
tion. These effects will also have to be taken into account 
when analysing upcoming Lyman-a experiments such as 
BOSS ( Slosar et al.||20iT ), which are expected to have sta- 
tistical errors an order of magnitude smaller than those of 
current data. 

There are many questions about high column density 
systems which remain unanswered by the present work. Al- 
though our AREPO simulations reduce the discrepancy be- 
tween DLA simulations and observations, completely elimi- 
nating it still seems to require stronger feedback processes. 
Furthermore, SPH simulations have historically failed to re- 
produce the velocity width of metal lines in DLAs, due to the 
suppression of mixing in SPH, which leaves metals concen- 
trated in small clumps, instead of being spread out through 
the halo ( |Tescari et al.|2009| ). Here the new AREPO code of- 
fers a lot of potential for progress. We intend to study these 
questions in future work, based on a more comprehensive 
model for feedback and metal enrichment (Vogelsberger et 
al., 2013, in preparation). 
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